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Common density-matrix functionals, the Muller and the power functional, have been benchmarked 
for the half-filled Hubbard dimer, which allows to model the bond dissociation problem and the 
transition from the weakly to the strongly correlated limit. Unbiased numerical calculations are 
combined with analytical results. Despite the well known successes of the Muller functional, the 
ground state is degenerate with a one-dimensional manifold of ferromagnetic solutions. The resulting 
infinite magnetic susceptibility indicates another qualitative flaw of the Muller functional. The 
derivative discontinuity with respect to particle number is not present indicating an incorrect metal¬ 
like behavior. The power functional actually favors the ferromagnetic state for weak interaction. 
Analogous to the Hartree-Fock approximation, the power functional undergoes a transition beyond 
a critical interaction strength, in this case however, to a non-collinear antiferromagnetic state. 
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I. INTRODUCTION 

Ab-initio calculations are dominated by density func¬ 
tional theory (DFT)^’^, which provides an efficient and 
accurate description of the electronic structure for most 
materials^. For materials with strong correlations, how¬ 
ever, many of the available density functionals yield poor 
results^’®. Most well known is the case of transition 
metal oxides, for which most density functionals pro¬ 
duce a qualitatively incorrect description®. However, 
also elementary chemical processes such as bond dissoci¬ 
ation are described poorly by currently available density 
functionals.® 

There is a quest to improve the description by bor¬ 
rowing from methods specifically designed for strongly 
correlated materials. Among them are LDA-|-U^, DFT- 
plus-dynamical mean-field theory*^^®, and DFT-plus- 
Gutzwiller approximation^ The guiding idea behind 
these approaches is to merge density functional theory 
with methods developed for the study of strong cor¬ 
relations for model Hamiltonians such as the Hubbard 
model. 

We consider reduced density-matrix functional theory 
(rDMFT)^^’^® to be a useful framework for a rigorous for¬ 
mulation of such hybrid theories.Reduced density- 
matrix functional theory can be viewed as a relative of 
DFT, which emphasizes orbital occupations rather than 
the density as basic variable. Such a description seems 
to be natural for correlated materials, because the latter 
are dominated by orbital physics. 

The link from rDMFT to many-particle wave func¬ 
tions has been established by Levy’s constrained-search 
algorithm^® on the one hand. The link to many-body 
perturbation theory and Green’s function, on the other 


hand, has been provided recently^® via the Luttinger- 
Ward functional^^. 

In order to avoid the full complexity of an explicit 
many-body description, most density-matrix function¬ 
als are not extracted from the exact expressions^®’^®. 
Rather, one proceeds analogously to the development of 
density functionals, namely by searching models^^^^® for 
the density-matrix functional, that capture the most es¬ 
sential physical effects while having an algebraic depen¬ 
dence on the density matrix. 

The development of such model density-matrix func¬ 
tionals relies on benchmark systems that allow one to 
evaluate their quality. Of particular interest are ex¬ 
actly solvable problems. Such studies have been per¬ 
formed for the Moshinsky atom^^, the homogeneous elec¬ 
tron gas^® and the Hubbard model^®’®®. Di Sabatino et 
al.®® performed an in-depth analysis of the method pro¬ 
posed by Sharma et al.®^ to evaluate the spectral function 
of the Hubbard dimer from the Muller density-matrix 
functional^®. 

As pointed out by Cohen, Sanchez and Yang®, many of 
the failures of current density functionals for correlated 
materials can be traced back to the derivative discontinu¬ 
ities present in a surprisingly simple system, namely the 
hydrogen or helium dimer in different charge states, i.e. 
HJ, H 2 , HeJ. Therefore, the two-site Hubbard model, 
the Hubbard dimer, can be considered as model system 
for the correlation effects present in a chemical bond. 

The most prominent failure of density functionals oc¬ 
curs during bond dissociation. If we denote the hopping 
parameter between the bonded atoms with t and the on¬ 
site interaction strength with U, bond dissociation is de¬ 
scribed by the limit t 0 at constant U. Thus, the 
system evolves from a weakly correlated state into the 
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strongly correlated limit U/t —>• oo as the bond is bro¬ 
ken. The large-interaction limit C/ —?► oo of the Hubbard 
model differs from the bond-dissociation limit only by 
the choice of the energy scale. 

One of the major arguments in favor of density-matrix 
functionals is that one of the most simple functionals, the 
Muller functionaH^, seems to provide a correct descrip¬ 
tion of the bond-dissociation problem, for which common 
density functionals fail^^. In this paper we study the 
performance of a class of commonly used model density- 
matrix functionals for the half-filled Hubbard dimer. We 
point out that, despite some successes, also these density- 
matrix functionals reproduce a number of features in 
a qualitatively incorrect manner. Thus, this work sets 
the stage for the development of entirely new class of 
functionals. 

In section H, we define our notation and introduce 
the basic concepts of density-matrix functionals. In sec¬ 
tion HI, we present the analytically exact treatment of 
the Hubbard dimer and describe its relevant properties. 
In section IV, we describe the numerical methodology 
of searching for the ground state for the model density- 
matrix functionals. In section V, we describe the results 
obtained with the Hartree-Fock approximation and the 
commonly used Muller and power functionals. In sec¬ 
tion VI we study the Hubbard dimer beyond the half 
filling and in section VH we discuss briefly the transfer- 
ability of our results to larger systems. Our results are 
summarized in section VHI. 


II. THEORETICAL FRAMEWORK 


A. General many-particle problem 


The many-particle Hamiltonian for interacting elec¬ 
trons can be expressed in terms of field operators '4’{x) 
and {x) in the form 

J d^X il^x) + Vextix)) 1 p{x) 


H = 


+ -z [ dPx [ d^x' •ip\x)'ip\x') -^- —'ip{x')'ip{x), 

2J J 4T:eo\r-r'\ 

( 1 ) 

where x = (r, a) is a combined position and spin variable. 
We use the shorthand f ddx = I integra¬ 

tion over positions and the sum over spin indices. The 
field operators obey the usual anticommutator relations 

tp^x),ip{x') = S{r- 
-I + 

A discrete, orthonormal one-particle basis set {Xa{x)} 
determines the creation and annihilation operators of 
electrons in the one-particle orbitals 


4 = y xAx)ip\x) 


Ca = J d^X X*a{x)4’{x). (2) 

In this one-particle basis set we obtain the discrete 


Hamiltonian 


Td = 2Z dapc, 


OL^ 


Co -b 




C^C ^C^C, 


t ■ 

aCp Cgv^ 


ol ^^5 


( 3 ) 


with the one-particle Hamiltonian 

= J d^x X*(x) + Vext{x)^ xp{x). (4) 

The off-diagonal elements of h are named hopping pa¬ 
rameters, and the diagonal elements are named orbital 
energies. 

The interaction matrix elements are 


t/< 




= d^ 


' 4 , c^Xa{x)x*p{x')x^{x)xs{x') 

47reo|f^— 


B. One-particle reduced density matrix 

The one-particle reduced density matrix of an ensem¬ 
ble of fermionic many-particle wave functions |$j) with 
probabilities Pj is defined as 

Pa/3 = 

j 

The density matrix is often represented by the eigen¬ 
values and eigenstates of the corresponding one-particle 
operator 

p = XI 

The eigenvalues of p are the occupations /„ and the eigen¬ 
states \(j)n) are named natural orbitals^^. Thus the den¬ 
sity matrix can be expressed by its eigenvalues and eigen¬ 
states in the form 

Pa./3 = '^{Xa\(t>n)fn{(l>n\Xp)- (8) 

n 

Not every hermitian matrix can be also be represented 
as the one-particle reduced density matrix of an ensem¬ 
ble of many-particle wave functions according to Eq. (6). 
A matrix that can be represented by an ensemble of 
fermionic N-particle wave functions is called ensemble N- 
representable. Coleman^^ has shown that eigenvalues of 
all ensemble N-representable one-particle reduced density 
matrices lie between zero and one and that all hermitian 
matrices with eigenvalues between zero and one are en¬ 
semble N-representable. 


C. Helmholtz potential and density-matrix 
functional 

The Helmholtz potential the thermodynamic 

potential for finite temperature and fixed particle num¬ 
ber, for a many-particle system can be expressed with the 
help of the density-matrix functional F'^ [p\ 
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Af}^N[h + W]= min statl ^ fniM^K) + Pf |(An)/«(</>« 

|<?^n),/ne 0,1 A,/i I 

^ n n 

M f ^ ^ fn '^1 ^ ^ ^m,n (^{4^n |(Am) ^m,n^ r 

\ n / m,n 


(9) 


where h = ^ |Xa)^a./3(X/3|- 

The reduce(i (density-matrix functional [p] is uni¬ 
versal in the sense that it depends only on the intrinsic 
properties of the electron gas, namely the interaction W, 
while it is independent of the one-particle Hamiltonian 
h. The chemical potential p is a, Lagrange multiplier 
that constrains the electron number to N. A^n are the 
Lagrange multipliers which enforce that natural orbitals 
I (Am) remain orthonormal. 

The reduced density-matrix functional 

Ff[p]=EH[p]+U,cAp] ( 10 ) 

is the sum of Hartree energy Eh and the exchange- 
correlation energy Uxc 

The Hartree energy Eh[p\ is obtained from the elec¬ 
tron density 


D. Hole function and the construction of 
density-matrix functionals 

1. Hole function 

In this section we discuss several exact properties of 
the hole function, which have been central to the develop¬ 
ment of density functionals, and in the following section 
we outline its role for the construction of approximate 
density-matrix functionals. 

The hole function h{r,r') allows to express the two- 
particle density in the form 

= n(f)n{r') + n{r)h{f,r'). (13) 

Note that the interaction-strength averaged hole function 
is used in DFT, while in rDMFT, the hole function at full 
interaction strength is of interest. 

The hole function integrates to minus one, 


(y a,^ 



(14) 


as 


Ph[p] 


1 

2 

1 

2 


^3^/ e^njrjnjP) 
47reo|r — r'\ 

^ ^ UcK^0^S,jP6,a.P'y^^- 



( 12 ) 


The exchange-correlation energy Uxc contains the com¬ 
plexity of the many particle problem. It is the elec¬ 
trostatic interaction of each electron with its exchange- 
correlation hole and the entropy term —TS. It should 
be noted that the exchange-correlation energy Exc of 
DFT also contains a contribution from the kinetic energy, 
which is absent in the quantity Uxc used in rDMFT. 

We restrict the present study to zero temperature and 
thus ignore the entropy term. To keep the notation sim¬ 
ple, we suppress the index for the inverse temperature in 
the remainder of the text. 

Having laid down the basic concepts and our notation, 
we proceed with the concept of the hole function as a 
tool for the construction of approximate density-matrix 
functionals. 


and it is always negative. 

These conditions constrain the shape of the hole func¬ 
tion strongly, so that the exchange-correlation energy can 
be predicted reasonably well already with simple assump¬ 
tions about the hole function. An insightful description 
of the hole function, which guided the development of a 
number of density-matrix functionals, has been given by 
Baerends and Buijse^®’^®. 

In the Hartree-Fock approximation, the hole function 
has the form 

nir) 

' ' m,n 

(15) 

As a consequence of the orthonormality of the natural 
orbitals, the sum rule Eq. (14) is obtained as 

/ dPr' h{pr') =flAuAAnix) = -1. (16) 

J ri{r) ^ 

The sum-rule is fulfilled exactly, when = /„ that is for 
integer occupations. For fractional occupations, however, 
the Hartree-Fock expression violates the sum-rule. 
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The exchange-correlation term in the Hartree-Fock ap¬ 
proximation is 

\p\ “ ~ 2 ^ / fmfn ^ ^ Uap,S'f 
m,n ot0^5 

><{Xl\(t>m){(t>m\Xa){X&\(l^n){(t>n\Xp)- (17) 


U, 


HF 


2. Construction of density-matrix functionals 


Most empirical density-matrix functionals maintain 
this general form of the Hartree-Fock exchange term, 

Uxc[p] — 2 ^ ^ ^ ^ 

m,n 

X (X7 I '('m) ('('m IXa) (X5 I (<)>n IX/?) , (18) 
but replace the factor fnfm in Eq. (17) by coefficients 
Cm,n with a different dependence on the occupations. 

Taking the hole-function in the Hartree-Fock approxi¬ 
mation Eq. (15) as a starting point, Mtiller^^ has shown 
that one can enforce the sum rule Eq. (14) also for frac¬ 
tional occupations with an ansatz 


/i(r,r') = 


1 


■ 4>*rn{x)4>n{x)(l)l{^)4>rn{x') 

m,n 


( 19 ) 

Muller identified p = 0 as the choice that minimizes the 
violation of the positive definiteness of the hole function. 
This is the value used in nearly all applications. 

Later, Sharma et al.^® invented the so-called power 
functional by introducing an additional parameter a in 
the dependence of the coefficients Cm,n on the occu¬ 
pations. They chose the form c^ „(a) = /„/“ that 
smoothly interpolates between the Miiller functional with 
a = \/2 and the Hartree-Fock approximation with a = 1. 
The main reason for this construction is according to 
Sharma et al.^® the well known overcorrelating behav¬ 
ior of the Muller functional, that will be mediated by a 
parameter a > 1/2. The coefficients Cm,n for the approx¬ 
imate density-matrix functionals considered in this work 
are summarized in table I. 


Hartree-Fock approximation 

HF — f f 
^m,n — JmJn 

Muller functionaF^ 

M _ f2 f2 
— JmJn 

power functional® 

(q'I = 

'-'m JmJn 


TABLE I. Dependence of the parameters Cm,n on the occupa¬ 
tions fn as defined in Eq. (18) for density-matrix functionals 
used in this work. 


III. HUBBARD DIMER 

The two-site Hubbard model, the Hubbard dimer, is 
the simplest model for the covalent bond and bond break¬ 


ing. 

The Hubbard dimer has a one-particle basis with four 
spin orbitals |xi,t); Ixi.l), IX2.t), |X2,i) , one for each site 
and spin. The only nonzero matrix elements of the one- 
particle Hamiltonian 

^a./3 = -til - ^R^,Ri3)5crc,a^ (20) 

are those with orbitals having the same spin Ua and up 
but different centers Ra and Rp. All nonzero elements 
have the value —t, where t is positive. The orbital ener¬ 
gies are chosen equal to zero. 

Also the interaction tensor has a simple form, namely 


Uap^jS 


U if a = 7 , /3 = (5 and Ra = Rp 
0 otherwise 


( 21 ) 


Thus, the Hamiltonian for the Hubbard dimer is 

^ = + (22) 

G 


with the interaction 




^i" ^ ^. 

-,G^i,G' ^i,G' ' 


(23) 


3. Total energy and density matrix 



U/t 


FIG. 1. Ground-state energy E of the half-filled Hub¬ 
bard dimer as function of interaction strength U/t for dif¬ 
ferent density-matrix functionals. Circles: Muller functional. 
Squares: power functional with a = 0.53. Triangles: Hartree- 
Fock approximation. Solid line: exact ground-state energy. 
The Muller functional produces the correct ground-state en¬ 
ergy at half filling. Non-magnetic states are indicated by open 
symbols and antiferromagnetic states by filled symbols. 

In Fig. 1, the total energy of the half-filled Hubbard 
dimer is shown as function of interaction strength, along¬ 
side with the results obtained from approximate density 
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matrix functionals. Some of these data have been pre¬ 
sented earlier^°. Here, we emphasize the ground states 
obtained without biasing the magnetic configuration. We 
follow the convention commonly adopted in the solid 
state community of showing the graph for varying in¬ 
teraction strength U and fixed hopping t, so that the 
hopping sets the energy scale. Considering the bond dis¬ 
sociation problem, the natural choice would be to keep 
the interaction strength constant, while reducing the hop¬ 
ping parameter. 

For the non-interacting case, i.e. at ?7 = 0, the wave 
function is a Slater determinant of bonding states with 
opposite spin 

= 0)) = i (cl^^ + 4,^) (4,^ + 4, J |0). (24) 

With \0) we denote the vacuum state. 

This wave function can be rewritten as superposition 
of two eigenstates of the interaction operator 

|$(17 = 0)) = - ^4,t^i.t + 

+ 2 (^14^2.4, ~ 1^)- 

The first wave function contains contributions with two 
electrons on the same site, i.e. ionic states. Its interac¬ 
tion eigenvalue is U. The second wave function describes 
two electrons with opposite spin on different sites. Its 
interaction eigenvalue is zero. 

The first term describes the double occupancy, that is 
the probability that two electrons are on the same site, 
which is penalized by the electron-electron interaction. 
The second term is attributed to left-right correlation, as 
it describes the probability that the two electrons are on 
different sites. 

As the interaction strength is increased, the contribu¬ 
tion of the first wave function, being responsible for dou¬ 
ble occupancy, is suppressed. The wave function obtains 
the form 

l^(^)) = + 4.t4, J 1^) cos + J) 

+ 1 ^) + 0 ■ 

(26) 

With a basis set in the order 
(IXi.t)Jxi 4 )JX 2 .t)JX 2 . 4 .)), the one-particle reduced 
density matrix has the form 

0 cos(2i?) 0 \ 

I 0 cos(2i?) 

0 10' 
cos(2i?) 0 1 / 

(27) 

The interaction energy is proportional to the double 
occupancy 

($(^?)|lT|$(^5)) = t/cos^ 0 (28) 


Pa,Pi'S) = 2 


1 

0 

cos(2i?) 

0 


and the non-interacting energy is 

($(i9)|h|$(z?)) =-2tcos(2i?). (29) 

The value of results from an equilibrium between 
the forces from the interaction energy Eq. (28) and those 
from the one-particle energy Eq. (29), which determines 
i?(17) as 


'd{U) = arctan 




TT 

4 ■ 


(30) 


The value 'i?(17) varies from zero to 7r/4 with increasing 
interaction strength. 

The resulting optimum energy has the form 



( 31 ) 


As the interaction increases, the wave function changes 
continuously from a Slater determinant of bonding states 
Eq. (25) at ?7 = 0 to a singlet state with antiferromag¬ 
netic correlations. During this process, the bond strength 
is weakened and the covalent bond vanishes completely 
in the limit of infinite interaction. This loss of covalent 
bonding can also be described as localization of electrons 
on opposite sites, which raises the kinetic energy as a 
consequence of Heisenberg’s uncertainty principle. 

What has been described here is what is called static 
correlation®: The states for finite interaction can no more 
be described by a single Slater determinant, but four 
Slater determinants are required. 


4- Natural orbitals and occupations 

Interestingly, the natural orbitals do not depend on 
the interaction strength U. They are the bonding and 
antibonding states 

l&>cr) := (IX 1 .<t) + 1x2 ,<t)) 

la,®") := (Ixi.<t) - 1X2, <t)) • (32) 

Both orbitals are spread over both atoms, and the natu¬ 
ral orbitals are identical to those of the non-interacting 
system. 

The loss of bonding is, however, expressed by the fact 
that the occupations become fractional. The occupations 
are shown in Fig. 2. Their exact values fb,(y for the bond¬ 
ing states and fa,a for the antibonding states are 

fb,a = \ + \ cos(2i9) 

/a,a = i - ^COs(2l?). (33) 

In the non-interacting case, the occupations are integer, 
with filled bonding states and unoccupied antibonding 
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states. In the limit of large interaction strength the oc¬ 
cupations approach i for all four natural orbitals. In 
this limiting case with equally occupied bonding and an- 
tibonding states, the net bond strength vanishes com¬ 
pletely. In the context of natural orbitals, we describe 
the effect as quantum fluctuations that create electron- 
hole pairs. These electron-hole pairs destroy the covalent 
bond with increasing interaction. 



FIG. 2. Occupations fb,cr and fa,a from Eq. (48) of the half- 
filled dimer obtained with the Muller functional as a function 
of interaction strength U/t. The striped regions indicate the 
range of occupations in the manifold of degenerate ground- 
states. The thick solid lines indicate the mean values for the 
pair of occupations in the corresponding striped region. It also 
represents the degenerate occupations for the non-magnetic 
solution of the Muller functional. The occupations of the 
non-magnetic solution of the Muller functional coincide with 
those of the exact ground state of the Hubbard dimer. 


5. Correlations 

In view of the following discussion, it is instructive to 
investigate the correlations of the electrons. The proba¬ 
bility for an electron to be on one site and the other on 
the other site, we name it “site correlation”, is given by 
the expectation value of 

C = (4,tC2.t - {4,1^2,1 - . (34) 

For a state where both electrons bunch on one site, the 
expectation value of this operator is one, while if they 
localize on opposite sites, the expectation value is minus 
one. If it is zero, then the electrons are delocalized, i.e 
there is no correlation between the positions of both elec¬ 
trons. The correlation operator C is a two-particle op¬ 
erator and is not accessible from the one-particle density 
matrix. The exact solution for the correlation expecta¬ 
tion value for the ground state is given by 

(C) =-sin(2i?(C/)), (35) 


where 'd{U) is given by Eq. (30). We can see in Fig. 3 that 
the site correlation vanishes without interaction, while 
the electrons anti-bnnch for strong correlation so that 
(C) approaches minus one. A site correlation of minus 
one indicates that each electron is fully localized either 
at one or at the other site, while the other is always 
at the other site. This is the basic notion of left-right 
correlation. 

Of interest will also be the magnetic nature of the wave 
functions. The operator for the spin on site i is 

- h ( ^ \ 

-S'* = 2 + *4.4.G.t • (36) 

V - 4,Ai j 

For the wave function in Eq. (26), the total spin ex¬ 
pectation value ((^i -I- iS' 2 )^) vanishes, and consequently 

the spin expectation value {Sf) on each site vanishes as 
well. Nevertheless, the spins on different sites are anti- 
ferromagnetically correlated, that is 

[1 +sin(2i?)] < 0. (37) 

o 

An antiferromagnetic correlation is already present in the 
non-interacting state, which expresses the non-vanishing 
contribution of the left-right correlated states to the 
Slater determinant built from bonding orbitals. As the 
interaction increases the left-right correlation doubles, 
which reflects in the increase of the antiferromagnetic 
correlation expressed in Eq. (37). 

IV. METHODS 

The natural orbitals and occnpations have been op¬ 
timized in the Car-Parrinello spirit^° using a fictitious 
Lagrangian of the form 


^ ^, 0-(y^nhfi,aag ,a F ( Qa,ra/(^ri)Qfl.n] 

n a,(3 n 



(38) 


The natural orbitals are given by the complex-valued co¬ 
efficients aa,n = iXal'Pn) aS 

I'Pn) — ^ ( \Xa)o-a,n (39) 

Ct 

and the occupations /„ = f{xn) are expressed by the 
real-valued dynamical variables Xn with f{x) = [1 — 
cos(a;)]/2. The orthonormality of the natural orbitals 
is enforced with the Lagrange multipliers A^.n, which 
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form a hermitian matrix, and the particle number is con¬ 
strained with the chemical potential 

In order to avoid any bias in our results, the wave func¬ 
tions and occupations are initialized as random numbers 
between zero and one. Then the constraints, i.e. or¬ 
thonormality of the natural orbitals and total particle 
number, are imposed. For the actual minimization, the 
Euler-Lagrange equations for the occupation variables Xn 
and the coefficients aa^n are propagated using the Verlet 
algorithm under the additional action of a friction term. 
The constraints are enforced with the help of Lagrange 
multipliers.^^ The friction term leads to energy dissipa¬ 
tion and monotonic decrease of the total energy until the 
ground state or a metastable state is reached. The phase 
space has been explored by repeating the calculation, in 
order to identify the global minimum and potential de¬ 
generate ground states. 

An analytical form of the natural orbitals has been ex¬ 
tracted by inspection of the results obtained numerically. 
The resulting ansatz for the natural orbitals has been ver¬ 
ified by optimizing the total energy in this subspace, and 
comparing the energies. While the numerical formulation 
is invariant under global spin rotations, spatial reflection, 
or application of a phase factor, the analytical results are 
given for a particular choice. 

V. PERFORMANCE OF DENSITY-MATRIX 
FUNCTIONALS 

A. Hartree-Fock approximation 


While the errors caused by non-spin polarized Hartree- 
Fock calculations are severe, they are not our main con¬ 
cern. Today’s electronic structure calculations should 
consider a spin-polarization whenever a magnetization 
provides a lower energy. Allowing for spin polariza¬ 
tion, i.e. as in unrestricted Hartree-Fock or spin-density 
functional theory, improves the description dramatically. 
Nevertheless, the transition from the weakly correlated to 
the strongly correlated regime still differs in many ways 
from the correct behavior. These differences are of inter¬ 
est in the following discussion. 

2. Antiferromagnetic solution 

If one allows for general variations of the density ma¬ 
trix, there is a crossover at U = 2t from a non-magnetic 
solution at small interactions to an antiferromagnetic so¬ 
lution at large interactions. 

One set of natural orbitals that describes the transition 
to the antiferromagnetic solution beyond U = 2t has the 
form 

1 ( 7 )) = + 1 t) cos( 7 ) -b I a, t) sin( 7 ) 
l^f^(7)) = +|64) cos( 7 ) - I a, 4 .) sin( 7 ) 
sin( 7 ) -b |a,t)cos( 7 ) 

= +l&4)sin(7) -b |a,4.)cos(7). (41) 

The hrst two natural orbitals are occupied and the re¬ 
maining two are unoccupied. 

The corresponding many-particle wave function. 


After having covered the main properties of the ex¬ 
act ground state of the half-filled Hubbard dimer in sec¬ 
tion HI, we turn now to the results obtained from ap¬ 
proximate density-matrix functionals. We begin with 
the Hartree-Fock approximation given in Eq. (17), which 
has been the starting point for the development of other 
density-matrix functionals investigated in this study as 
discussed in Sec. HD2. 


1. Non-magnetic solution 

If one constrains the density matrix to remain non¬ 
magnetic, the natural orbitals do not depend on the in¬ 
teraction strength. The corresponding total energy has 
the form 

E^^{U) = -2t+^. (40) 

The energy Eq. (40) for the Hubbard dimer with two 
infinitely separated atoms, that is in the limit of t —>■ 0 , 
results in a non-zero energy while the correct en¬ 
ergy vanishes, because each isolated atom has a single 
electron that does not interact with itself. This reflects 
the well known difficulty of restricted, i.e. non-spin- 
polarized, Hartree-Fock to describe the dissociation of 
chemical bonds. 


|ci>^^(7)) = 


X 



-b C2 I COS 

-4.iCOS 


(^fo): 

( 7 - 4 )] |o), 


(42) 


is a single Slater-determinant in the basis of the natural 
orbitals. For 7 = 0, we recover the ground state of the 
non-interacting limit given in Eq. (24). 

The many-particle wave-function Eq. (42) has the 
one-particle reduced density matrix in the basis 

(IXi.t)>IXi.4)JX2.T)jA2.i)) 

l-bsin( 27 ) 0 cos( 27 ) 0 \ 

0 1 —sin(27) 0 cos(27) 

cos( 27) 0 I —sin(27) 0 

0 cos( 27 ) 0 l-bsin( 27 ) J 

(43) 

The interaction energy is, 

($^^(7)|IF|4>^^(7)) = ic/cos2(27) (44) 

and the non-interacting energy is given by 

(<j)^^( 7 )|/i|(j)-f^-f’( 7 )) = —2tcos(27). (45) 









Increasing the parameter 7 in the wave function from 0, 
i.e. the non-interacting limit, allows one to trade part of 
the covalent bond, i.e. the kinetic energy, for a reduction 
of the interaction energy. 

The total energy is minimized by 


l{U) 


0 for C7 < 2t 

^ arccos (|^) for U > 2t. 


(46) 


For U < 2t, the system remains non-magnetic and the 
natural orbitals are given by bonding and antibonding 
orbitals as in the case of non-magnetic dimer. But for 
U > 2t, the system becomes an antiferromagnet, whereas 
the exact many-particle wave function is a singlet with 
antiferromagnetic correlations. The antiferromagnetic 
state is a superposition of a singlet and a triplet wave 

function and thus it is not an eigenstate of S'^. We can 
paraphrase it as a violation of rotational symmetry in the 
spin degrees of freedom, i.e. of SU(2) spin symmetry. 


B. Muller functional 


the bonding states are occupied, the remaining terms of 
Eq. (47) vanish and the Muller functional gives the same 
result as the Hartree-Fock approximation. 

The occupations are obtained as the minimum of 
Eq. (47) for occupations between zero and one that add 
up to the total particle number of fV = 2. For a given 
interaction strength, we find that the minimum condition 
does not define a single point, but a line of degenerate 
states parameterized by the parameter s 


/a“(5) 


1 


1-f i?2 

1 

1-f 


= 


1 


1-f i?2 
1 

1 -bi?2 



(48) 


where R = At/U + \Jl + (4t/{7)^. 

The requirement, that the occupations remain between 
zero and one, limits the parameter s to the interval 


Muller’s approximation to the density-matrix func¬ 
tional introduced in Sec. IID 2 leads to the exact ground- 
state energy for the half-filled Hubbard dimer for all in¬ 
teraction strengths^®’^^. In contrast to the Hartree-Fock 
approximation, there is no unphysical transition to an 
antiferromagnetic state. 


1. Magnetic solutions 


[ R^il + R'^)' R'^il + R^)\' ^ ' 

The range of the occupations, which minimize the total 
energy Eq. (47), is shown in Fig. 2 as a function of inter¬ 
action strength Ujt. In the limit of infinite interaction 
strengths, we have i? = I respectively s € [—1/2,1/2] 
and the possible occupations ^(s) = 1/2-\- as cover 
the whole range from zero to one. 


Even though the Muller functional produces exact 
ground-state energies for the half-filled Hubbard dimer, 
we also detected a major flaw, namely that there is a one- 
dimensional manifold of magnetic states which are degen¬ 
erate to the correct non-magnetic solution. The infinite 
magnetic susceptibility obtained with the Muller func¬ 
tional is in contrast to the exact behavior: At zero tem¬ 
perature and finite interaction strength, the true mag¬ 
netic susceptibility vanishes because of the finite singlet- 
triplet splitting.'‘ 3 d 4 

Our unbiased optimizations using the Muller func¬ 
tional result in natural orbitals equivalent to the exact 
ones given in Eq. (32), namely the bonding and anti¬ 
bonding orbitals. 

With the natural orbitals of Eq. (32), the total energy 
for the half-filled dimer obtained from the Muller func¬ 
tional can be expressed solely by the occupations as 

(47) 

The first two terms, which are independent of the oc¬ 
cupations, are identical to the total energy Eq. (40) of 
the spin-restricted Hartree-Fock approximation. If only 



FIG. 3. Site correlation {C) as defined in Eq. (34) of the 
half-filled Hubbard dimer as a function olU/t. With increas¬ 
ing interaction strength U/t the site correlation shows the 
transition from delocalized electrons (C) = 0 to the left-right 
correlated state with (C) = — 1 . 

All solutions, except the physical one with equally oc¬ 
cupied bonding states and equally occupied antibonding 
states, have a magnetic moment. Hence, the magnetic 
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susceptibility predicted by the Muller functional is infi¬ 
nite for all finite interaction strengths. 

The magnetization of each site in the ground state 
of the Miiller functional obtained with the occupations 
given by Eq. (48) has the form 

m^s) = i [/b^(s) - /fc^(s) -b /i^(s) - /a*l(s)] 

= {l + R^)s^lB (50) 

with the Bohr magneton ^b- It can assume any value 
with \mz\ < 1/i?^ ^iB■ 

In Fig. 4, the density-matrix functional of Muller 
is compared to the exact density-matrix functional in 
the range of degenerate ground states of the Muller 
functional. The exact density-matrix functional is ob¬ 
tained from a constrained search over an ensemble of 
fermionic many-particle wave functions^® for density ma¬ 
trices parametrized by Eq. (32) and Eq. (48). The enor¬ 
mous difference in the functionals illustrates the severe 
problems of the Muller functional to describe the mag¬ 
netic structure properly and indicates a systematic flaw 
in the functional. 



FIG. 4. Muller density-matrix functional (dashed line) and 
the exact functional (solid line) as a function of the line pa¬ 
rameter s for U = 4t. The density matrix p(s) is given by 
Eq. (32) and Eq. (48). The values of the exact functional 
have been obtained by a constrained search over an ensem¬ 
ble of many-particle wave functions. The point s = 0, where 
the Muller functional and the exact functional coincide, cor¬ 
responds to the symmetric solution (niz = 0). 


2. Ojf-site interaction 


also an electron-electron interaction V between the sites 
W = - Ucj ^rCi^a-'Ci^cr 


9 X/ X/ ■ 


(51) 




Using the density matrices from the degenerate man¬ 
ifold of ground states without offsite interaction, i.e. 
with bonding and antibonding states as natural orbitals 
Eq. (32) and the occupations from Eq. (48), the effect of 
the off-site interaction has been explored up to first order 
in the off-site interaction V. This leads to 

E^[V,s\ = E^[0,s] + 


X 


2 

l-bi?2 



+ 0(V^). 


(52) 


E^[0, s] is the s-independent total energy obtained with 
the Muller functional for the Hubbard dimer in the ab¬ 
sence of an offsite interaction. It is given by Eq. (47) and 
Eq. (48). 

As shown in Fig. 5, the off-site term lifts the degener¬ 
acy of the ground states of the Muller functional. The 
non-magnetic solution with s = 0 is now favored. This 
indicates that this artificial degeneracy may not be im¬ 
mediately apparent in real systems. 

Nevertheless, as evident from the comparison with the 
exact functional shown in Fig. (4), the changes produced 
by the off-site term are far too small: In order to pro¬ 
duce an energy difference between the maximally polar¬ 
ized state (see Eq. 49) and the unpolarized state compa¬ 
rable to the exact result shown in Fig. 4, an unrealisti¬ 
cally large offsite interaction parameter of order V = lOt 
would be required. 



- 0.02 - 0.01 0 0.01 0.02 


line parameter s 


The finding of an infinite magnetic susceptibility raises 
the question whether this finding transfers to more real¬ 
istic systems. One of the major restrictions of the Hub¬ 
bard model is the limitation to pure on-site interactions. 
Therefore, we extended the Hubbard model to include 


FIG. 5. Energy AA = E[V, s] - E[V, s = 0] of Eq. (52) of the 
Hubbard dimer obtained with the Muller functional including 
an off-site interaction in first-order perturbation theory with 
U — 4:t along the manifold Eq. (48) of ground states. The 
point s = 0 indicates the non-magnetic solution. 
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C. Power functional 

After having investigated the Hartree-Fock approxi¬ 
mation and the Muller functional, we consider now the 
power functional invented by Sharma et al.^®, which we 
described in Sec. IID 2. 



U/t 


FIG. 6. Occupations fi as a function of U/t for the power 
functional with a = 0.53. Solid dots have been obtained from 
an unbiased optimization of the power functional. The solid 
lines are obtained from a restricted optimization using the 
non-collinear natural orbitals of the ansatz Eq. (56). The 
diagonal crosses are obtained from a restricted optimization 
using the collinear natural orbitals Eq. (55) of the Hartree- 
Fock approximation. 


The occupations of the Hubbard dimer obtained from 
the power functional are shown in Fig. 6 for a value 
a = 0.53. Whereas the density-matrix functional in the 
Hartree-Fock approximation produces integer and pair¬ 
wise identical occupations, the power functional produces 
fractional occupations which are not identical in pairs. 

Near U = 6t, we observe a transition. This transition 
separates the Muller-like behavior at small interactions 
from a Hartree-Fock-like behavior at large interactions. 

• At small interactions, the solutions are analogous 
to those of the Muller functional. However, from 
the manifold of degenerate ground states of the 
Muller functional, the power functional favors the 
state with maximal ferromagnetic moments. 

• At larger interactions, the ground state under¬ 
goes a transition into a non-collinear ground state. 
For very large interaction the state approaches the 
Hartree-Fock-like antiferromagnetic state. 


exists a manifold of degenerate ground-state density ma¬ 
trices on the line given by Eq. (48). If we increase the 
parameter a of the power functional infinitesimally as 
a = ^ + e where e > 0, and restrict ourselves to interac¬ 
tion strengths U/t where the natural orbitals are bonding 
and antibonding states, Eq. (32), the total energy along 
the line given by Eq. (48) is 


<=i+,(5) = 2t 


U 


1 


-l]+U 


-tE 


<7—:kl 


1 


(53) 

l+2e 




+ as 


where R = At/U + ^/T+J^t/U)^. The energy in Eq. (53), 
shown in Fig. 7, has a negative curvature along the line 
parameter s and the minima lie at the boundaries given 
in Eq. (49). 

At these boundaries, the extreme non-symmetric solu¬ 
tions of the Miiller functional, one of the states is always 
fully occupied (See Fig. 2) because this maximum oc¬ 
cupation limits the range of degenerate solutions. This 
explains the corresponding observation in Fig. 6. 

Unfortunately, any change of the parameter a away 
from the value of the Miiller functional, destroys the non¬ 
magnetic ground state in favor of an unphysical ferromag¬ 
netic state. 



line parameter s 


FIG. 7. Total energy difference AE = Ra=i/ 2 +e [P(s)] - 
E^^i/ 2 [p(s)] given by Eq. (53) for U = At using the power 
functional approximation with e = 10~^ as a function of the 
line parameter s that parametrizes the one-particle reduced 
density matrix according to Eq. (48). 


2. Large-interaction regime 


1. Ferromagnetic solution in the weakly interacting regime: 

The occupations in the weakly interacting regime can 
be understood as follows: In case of the Muller density- 
matrix functional, we have shown in Sec. VB that there 


The Hartree-Fock approximation exhibits a transition 
from a non-magnetic state to an antiferromagnetic state 
at U = 2t. This transition is absent in the Muller func¬ 
tional, but it is present in the power functional for all 
other values of a > 
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In order to explore, how the power functional inter¬ 
polates between these two extreme cases, we calculated 

the product (Si) ■ (S 2 ) of the spin expectation values at 

the two sites of the dimer. A positive value of (^i) • {S 2 ) 
corresponds to a ferromagnetic, a negative value to an 
antiferromagnetic spin alignment. The maximum abso¬ 
lute value is /4. 

For the Hubbard dimer at half filling, (^i) • {S 2 ) is 
shown in Fig. 8 as function of interaction strength U/t 
and the parameter a of the power functional. For the 
Muller functional discussed in Sec. VBl, , i.e. for 
a = 1/2, we consider the solution with the strongest 
polarization, because this is the state that continuously 
matches to the solutions of the power functional. In this 

ferromagnetic state, (Si)-{ 82 ) is positive. Unfortunately, 
the correct non-magnetic state is not a ground state of 
the power functional for a > i. 

At a critical interaction strength Uc{a) the power func¬ 
tional exhibits a transition from this ferromagnetic state 
into a complex non-collinear state with a mostly anti¬ 
ferromagnetic spin alignment. The angle between the 
magnetization on the two sites is shown in Fig. 9. 

Fig. 8 clearly shows the location of the transition be¬ 
tween the ferromagnetic and the antiferromagnetic non- 
collinear regime. The critical interaction strength 14 (a) 
of this transition is infinite for the Muller functional. 
As the parameter a is increased, the critical interac¬ 


tion strength falls off rapidly and approaches the value 
?7c(a = I) = 2t of the Hartree-Fock approximation. 

Thus, the power functional exhibits a Hartree-Fock- 
like transition into an antiferromagnetic ground state ex¬ 
cept for the limiting case, the Muller functional. By 
choosing the parameter a sufficiently close to 1 / 2 , the 
transition can be shifted into a regime that is physically 
less important. 

a. Collinear approximation using the Hartree-Fock 
natural orbitals In order to get a qualitative under¬ 
standing of the asymmetric occupations (Ref. Fig. 6 ) 
and the critical value of interaction strength Uc of the 
transition to antiferromagnetic solutions (Ref. Fig. 8 ), 
we use an ansatz that covers both extreme cases, namely 
the Muller functional with a = ^ and the Hartree-Fock 
approximation with a = 1. These are, one the one hand, 
the asymmetric natural orbitals Eq. (41) that can de¬ 
scribe the antiferromagnetic state of the Hartree-Fock 
approximation. On the other hand, the ansatz allows 
for fractional occupations to capture the nature of the 
ground state of the Muller functional. 

With this ansatz, the one-particle reduced density ma¬ 
trix p(/i,..., / 4 , 7 ) is a function of occupations /„ and 
the angle 7 and the corresponding total energy E^’°‘ ob¬ 
tained with the power functional is 

EaHfu ■ ■ ■) /4, 7 )] = 7 )] 

+ Ff[p(/i,...,/4,7)] (54) 

where 


E = -tcos( 27 ) (/i -f /2 - /a - fi) 

F^Pifu- ■ •, /4, 7 )] = ^ [(/i + /2 + /3 + Uf - (/r + - (/2“ + 

U 


+ J sm^(2j) [if 2 + /3 - /4 - hr - ifl - f^r - {f 2 - Ia? 


(55) 


An approximation, which is a strict upper bound, for 
the total energy of the power functional is obtained by 
minimizing Eq. (55) for a half-filled system with occupa¬ 
tions between zero and one. 

As a characteristic example, the resulting occupations 
for a = 0.53 are shown in Fig. 6. The properties of this 
ansatz with regard to the description of the transition 
to the antiferromagnetic state will be investigated in the 
following section after a more general discussion of the 
transition. 

The ansatz using the collinear natural orbitals Eq. (41) 
and arbitrary occupations is, however, not able to de¬ 
scribe the true ground state for the power functional in 
the strongly interacting regime. The energy difference of 
the ansatz to the unbiased solution is shown in Fig. 10. 
The deviation is largest near the transition. The tran¬ 
sition point is slightly displaced by the collinear ansatz, 
which explains the sharp rise. For larger interactions, the 


error due to the collinear approximation falls off rapidly. 
It should be noted that the overall error due to the re¬ 
stricted ansatz is apparently small. For the parameter 
a = 0.53 used in Eq. (10), the maximum error in the 
energy is less than 1 % of the binding energy. 


b. Beyond the collinear approximation The ansatz 
using the Hartree-Fock natural orbitals already gives a 
fairly good description of the ground state of the power 
functional. How do the natural orbitals of the power 
functional differ from those of the Hartree-Fock solution? 


In the large interaction region, the power functional 
produces non-collinear ground states. The natural or¬ 
bitals of the power functional can be represented as su- 
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FIG. 8. Scalar product {Si) ■ (S 2 ) of the spin expectation vec¬ 
tors on the two sites of the Hubbard dimer as an indicator for 
the transition to the antiferromagnetic state within the power 
functional approximation with the parameter a for the Hub¬ 
bard dimer at various interaction strengths. A positive value 
indicates a ferromagnetic state, a negative value an antiferro¬ 
magnetic state. For the Muller functional, i.e. a = 1/2, the 
dashed line represents the result for the symmetric solutions 
and the solid line the corresponding degenerate result for the 
degenerate maximally polarized state. 



u/t 


FIG. 9. Angle between the spin expectation vectors (Si) and 
(S 2 ) on the two sites of the Hubbard dimer as function of the 
interaction strength U. Dashed line: power functional with 
the parameter a = 0.53; solid line: Hartree-Fock approxi¬ 
mation. This angle is a measure of collinearity of natural 
orbitals. 

perpositions of bonding and antibonding states, 
l<^f) = l^>t) cos(^i) - |a,4.) sin(/3i) 

l<^^) = \b,i) cos(^2) - |a,t) sm{^2) 

l</>r) = t) sin(/3i) -h |a, i) cos(/3i) 

l</>4) = 1^, i) sin(^2) + |a, t) cos(/32) ■ (56) 

The two angles /3i and ,02 are free variational parame¬ 
ters. The natural orbitals of the non-interacting system, 
respectively those of the Muller functional are obtained 


with j3i = (32 = 0. The values of the two parameters are 
shown in Fig. 11 for one example of the power functional. 

In the Hartree-Fock approximation, respectively in the 
power functional with the collinear ansatz, the pair of 
bonding and antibonding orbitals that contribute to a 
natural orbital, given in Eq. (41), have the same spin di¬ 
rection. This results in the localization of the electron on 
one or the other site of the dimer. The emerging picture 
is appealing because it reflects the left-right correlations 
of the electrons. The admixture of antibonding states 
to the natural orbitals for the two spin directions is the 
same. Thus, there is no symmetry-breaking charge dis¬ 
proportionation. 

The natural orbitals Eq. (56) of the power functional 
are composed of bonding and antibonding orbitals with 
opposite spin directions. This leads to natural orbitals 
with equal weight on both sites, but the spins on both 
sides have a finite angle between them. The state has an 
intrinsically non-collinear, even though still a coplanar 
spin structure. 

The admixture of antibonding states in the two pairs is 
independent in the power functional, so that the natural 
orbitals contain two independent free parameters, namely 
(3i and (32- 

The net magnetic moment of each of the four natural 
orbitals points along the same direction. For the choice 
in Eq. (56), this is the z-direction. The parameters (3i 
and /?2 control the relative angles of the local moments on 
the two sites of the dimer for each of the natural orbitals. 
This angle is 4/3i for the orbitals |^f) and \(j)^) and it is 
4/?2 for the orbitals |0^) and \4>^). The natural orbitals 
are pairwise antiparallel: On any given site \4>i) and \((i^) 
have local moments in opposite directions. Similarly, this 
holds for \(j) 2 ) and \4>^). 

It seems that the ground states of the power functional 
do not connect continuously to those of the Hartree-Fock 
approximation, because the natural orbitals belong to 
different classes. This is, however, not so: The ansatz for 
the natural orbital Eq. (56) connects smoothly to those 
of the Hartree-Fock approximation in Eq. (41) when the 
two parameters (3i and (32 become equal, and further¬ 
more the occupations become integer. This limit of the 
ansatz Eq. (56) for the power functional describes, how¬ 
ever, an antiferromagnet with the local moments aligned 
along the x-direction, while the ansatz of Eq. (41) for the 
Hartree-Fock solution is polarized along the z-direction. 
Thus they are related by a global spin rotation, which is 
a symmetry of the Hamiltonian. 


VI. BEYOND HALF-FILLING 

Up to now, we considered only the half filled case of the 
Hubbard dimer. Here we consider also deviations from 
the particle number N = 2. 

To avoid mathematical complications, we define E{N) 
thermodynamically consistent as the zero-temperature 
limit of the Helmholtz potential /3 —^ 00 , which in turn 
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FIG. 10. Energy difference AE of the power functional with 
a — 0.53 between the density matrices obtained by a con¬ 
strained and an unbiased optimization. AE for the case 
of constrained optimization with the natural orbitals of a 
Hartree-Fock in Eq. (41) is the solid line with circle symbols, 
while the AE obtained from constrained optimization with 
the non-collinear natural orbitals of Eq. (56) is the dashed 
line with square symbols. 



U/t 


FIG. 11. Parameters Pi and P 2 defining the natural orbitals 
Eq. (56) of the power functional for a = 0.53 as function of 
the interaction strength. 

is constructed from the grand potential by a Legendre- 
Fenchel transform 

E{N) = lim max 

^—¥OQ fJ. 

(57) 

The trace is performed over the fermionic Fock space. 

It can be shown that the total energy E{N) consists of 
piecewise linear segments between integer particle num¬ 
bers. Thus the slope of the total energy E{N), the chem¬ 
ical potential fi = dE/dN, is usually^® discontinuous at 
integer occupations. This derivative discontinuity gives 
the fundamental band gap which is defined as the dif¬ 
ference of electron affinity and ionization potential. The 


band gap is relevant, not only as an estimation related to 
optical spectra, but, more importantly, for the response 
functions and chemical equilibria. Therefore, we investi¬ 
gate whether the derivative discontinuities are properly 
described by the approximate density-matrix functionals. 

The total energy E{N) of the exact solution and sev¬ 
eral power functionals is shown in Fig. 12 and the corre¬ 
sponding chemical potential in Fig. 13. For the Hubbard 
dimer, the derivative discontinuity at = 2 is due to a 
combination of the one-particle gap and the interaction. 
The derivative discontinuity at iV = 1 is, however, en¬ 
tirely due to the interaction. These features are clearly 
visible for the exact calculation shown in Fig. 12. 

In the Hartree-Fock approximation, the energy for 
fractional occupations has a negative curvature for 1 < 
N < 3. As a result, the derivative discontinuities are 
larger than in the exact solution. It reflects the well 
known observation that Hartree-Fock calculations over¬ 
estimate band gaps. This observation can be rationalized 
with a lack of screening in the Hartree-Fock approxima¬ 
tion that reduces the effective interaction strength. 

The Muller functional, however, fails to give any 
derivative discontinuity. It is surprising, that a functional 
that performs as well as the Muller functional for iV = 2 
is completely unable to capture the correct physics be¬ 
yond half filling. It adds to the simplified picture that 
the Muller functional behaves very metal-like: It does 
not have a band gap and and its magnetic susceptibility 
is infinite. 

Except for the Hartree-Fock limit, also the power func¬ 
tional lacks a derivative discontinuity. This is apparent 
from Fig. 13. For small a, that is the Miiller-like regime, 
the power functional behaves analogous to the Muller 
functional itself. In the parameter regime of the antifer¬ 
romagnetic ground state, however, the chemical potential 
makes a continuous transition between two distinct linear 
functions of the particle number. 

This behavior of the power functional for the Hubbard 
dimer is analogous to that observed earlier for finite'*’®’^^ 
and extended systems^®’^®. 

In order to extract values for the band gap despite 
of the absence of a derivative discontinuity, Sharma et 
al.^® proposed the extrapolation method, which exploits 
the behavior of E(N) further away from the Fermi level. 
Sharma et al. exploit that the chemical potential makes 
a transition between two linear functions. The extrap¬ 
olation of these linear functions to the integer particle 
number yields an offset which is identified with the band 
gap. This method yields finite band gaps in the appropri¬ 
ate parameter range of the power functional, where the 
Muller functional incorrectly predicts a vanishing band 
gap^®. Surprisingly, the band gaps obtained using the ex¬ 
trapolation method from the power functional agree well 
with experimental results even for non-magnetic calcula¬ 
tions. 

Our results for the Hubbard dimer shown in Fig. 13 
indicate that the band gap obtained with the extrap¬ 
olation method^® can be tuned between zero and the 


In (ti- 
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Hartree-Fock value by adjusting a. Signatures of this 
behavior have been observed in studies that investigated 
the dependence on the power functional parameter a for 
realistic systems^®’^^. 

The absence of a true derivative discontinuity using the 
power functional and the tunability of the band gap de¬ 
termined with the extrapolation method is not limited to 
the antiferromagnetic ground state. As shown in Fig. 14, 
the Hubbard dimer behaves qualitatively similar, when 
the spin polarization is suppressed. In the nonmagnetic 
calculations, the onset of a finite band gap obtained by 
the extrapolation method is delayed to larger power pa¬ 
rameters a. This finding is analogous to that observed for 
NiO, for which non-magnetic calculations find a metal¬ 
lic ground state for a < 0.65^®, whereas non-collinear 
calculations find an insulating ground state already for 
a = 0.56^°. 



FIG. 12. Ground-state energy E{N) of the Hubbard dimer 
with U = 5t in units of the hopping parameter t as function of 
particle number N. The critical power functional parameter 
for the transition to an antiferromagnetic state lies at a « 
0.54 for the given interaction strength. Dashed line: exact 
solution, crosses: Hartree-Fock approximation, open circles: 
power functional with a = 0.7, filled circles: power functional 
with a = 0.58, triangles: power functional with a = 0.53, 
squares: Muller functional. 


VII. BEYOND THE DIMER 

The question remains whether the findings for the 
Hubbard dimer persist in larger systems with more de¬ 
grees of freedom. This is relevant for calculations of more 
complex systems having large unit cells. For this purpose 
we performed calculations for the power functional for 
Hubbard rings and Hubbard chains. 

Figure 15 shows the occupation numbers for a half- 
filled Hubbard ring at an intermediate interaction 
strength of C/ = 5t, which like the Hubbard dimer, has 
an antiferromagnetic ground state in the Hartree-Fock 



N 

FIG. 13. Ghemical potential ^{N) of the Hubbard dimer with 
17 = 5t in units of the hopping parameter t as function of 
particle number N. The behavior of the power functional 
with 1/2 < a < 1 close to half-filling is shown in the inset. 
Dashed line: exact solution, crosses: Hartree-Fock approxi¬ 
mation, open circles: power functional with a = 0.7, filled 
circles: power functional with a = 0.58, triangles: power 
functional with a = 0.53, squares: Muller functional. 



FIG. 14. Ghemical potential ^{N) of the Hubbard dimer with 
17 = 5t in units of the hopping parameter t as function of 
particle number N, when the density matrix is restricted to 
be non-magnetic. The behavior of the power functional with 
1/2 < a < 1 close to half-filling is shown in the inset. Dashed 
line: exact solution, crosses: Hartree-Fock approximation, 
open circles: power functional with a = 0.95, filled circles: 
power functional with a = 0.9, triangles: power functional 
with a = 0.85, squares: Muller functional. 


approximation. For the Miiller functional we obtain frac¬ 
tional occupations as for the dimer. While the fractional 
occupations deviate from the exact result, their deviation 
from integer occupations are of the same order of mag¬ 
nitude as in the exact solution. The power functional 
exhibits abrupt transitions to an antiferromagnetic state 
around Oc ~ 0.58 very analogous to the Hubbard dimer. 
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FIG. 15. Occupations of the half-filled six-site Hubbard-ring 
with U = 5t for the power functional as function of the param¬ 
eter a (solid lines). The dashed horizontal lines indicate the 
occupations of the exact many-electron description. Evident 
are the rather abrupt transitions from fractional to integer 
occupations. 



FIG. 16. Occupations of the six-site Hubbard-chain with 
seven electrons and U = 5t for ground states of the power 
functional as function of the parameter a spanning the range 
from the Muller functional (q = 1/2) to the Hartree-Fock 
approximation (a = 1). Evident are the rather abrupt tran¬ 
sitions from fractional to integer occupations. 


For a six-site Hubbard chain with seven electrons, i.e. 
one electron more than half filling, the pattern of tran¬ 
sitions is even more complex: This behavior is shown in 
Figure 16. There are now three transitions: 

1. A continuous transition between a ~ 0.567 and 
a ~ 0.569 from the non-magnetic Muller ground 
state to a state with collinear spins in the pattern 
tittitj which is only stable in a small window of 
parameters. 

2. Around a ~ 0.576 there is a non-smooth transition 
to a state with an antiferromagnetic spin-structure. 


i.e. tJ-tJ-tJ- 

3. Beyond a ~ 0.75, the antiferromagnetic structure 
breaks up and evolves into the HF-ground state 
having a spin-structure given by tittJ-t- 

These examples demonstrate that the power functional 
can generate a variety of magnetic states even for simple 
systems. 


VIII. CONCLUSION 


The popular density-matrix functionals, the Muller 
functional^^, the Hartree-Fock approximation and the 
power functionaP^, which continuously interpolates be¬ 
tween the other two, have been benchmarked for the 
Hubbard dimer. The Hartree-Fock approximation is, 
for the Hubbard rnodel^^^^®, analogous to hybrid den¬ 
sity functionals®^, that admix a portion of exact ex¬ 
change to the exchange-correlation energy. The local in¬ 
teraction of the Hubbard model acts analogous to the 
range separation®^’®®, which suppresses the long-ranged 
Coulomb interaction in the Fock term. In this respect, 
the Hartree-Fock approximation also captures the main 
effects of the LDA-bU method^. 

Particular emphasis has been given to left-right corre¬ 
lation, the dominant correlation effect for bond dissocia¬ 
tion, which is not captured in local density functionals®. 
Left-right correlation describes that electrons localize on 
opposite sites of the dimer. This electron correlation, 
which increases with the interaction strength, avoids the 
energetic cost of the Coulomb repulsion due to double 
occupancy of a site. In the Hartree-Fock approximation, 
this left-right correlation leads to an antiferromagnetic 
state with a spin-up electron mostly localized on one side 
and the spin-down electron on the other. This so-called 
broken-symmetry state disagrees with the exact solution, 
which is a singlet state, having no local moments, but 
nevertheless antiferromagnetic correlations similar to the 
broken symmetry state. More importantly, however, the 
antiferromagnetic transition is an abrupt one and not a 
continuous buildup of antiferromagnetic correlations as 
in the exact solution. The result is a qualitatively incor¬ 
rect shape of the total energy during bond dissociation. 

The Muller functionaP^ establishes left-right correla¬ 
tion in a fundamentally different manner: while the nat¬ 
ural orbitals are mostly - in the Hubbard dimer exactly- 
independent of the interaction, the occupations become 
fractional, which reflects the creation of electron-hole 
pairs that screen the interaction. One of the main suc¬ 
cesses of the Muller functional besides being able to pro¬ 
duce fractional occupations correctly, is that it captures 
the continuous nature of the transition to the left-right- 
correlated state. 

Our calculations avoid any bias and allow for arbitrary 
non-collinear spin-polarized states. This strategy shall 
bring all potential problems to the surface, that would 
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be present in large scale electronic structure calculations 
using these density-matrix functionals. 

Our first observation is that the ground state for the 
Muller functional, which does not have local moments, is 
degenerate with a one dimensional manifold of ferromag¬ 
netic states. Thus the dimer has infinite magnetic sus¬ 
ceptibility when described with the Muller functional, in 
contrast to the vanishing zero-temperature susceptibility 
of the exact solution of the Hubbard dimer. This large 
magnetic polarizability is likely to cause severe problems 
in extended electronic structure calculations. 

When turning to the power functional^®, we find that 
the system behaves analogous to the Muller functional 
for small interactions, while it exhibits a transition to a 
Hartree-Fock-like antiferromagnetic state for large inter¬ 
actions. The critical interaction, where this transition oc¬ 
curs, drops rapidly with increasing a from infinity in the 
Muller functional to the Hartree-Fock value Ucrit = 2t. 

In the small-interaction regime the system is weakly 
pinned in the ferromagnetic state corresponding to the 
largest moment of the ground-state manifold of the 
Muller functional. 

Our calculations indicate a major deficiency in the de¬ 
scription of magnetic properties for this class of density- 
matrix functionals. The problems persist in modified 
form also for more general Hamiltonians, which include 
off-site Coulomb interactions, and for more extended sys¬ 
tems. 

Besides the bond-dissociation problem, we investigated 
the derivative discontinuity®’®^ with changing the number 
of electrons. A balanced description of the electron affin¬ 
ity and ionization potential is essential for a qualitatively 
correct description of charge transfer. We find that the 
metal-like behavior of the Muller functional persists: The 
discontinuity of the exchange-correlation energy even off¬ 
sets the one of the kinetic energy. The Muller functional 
describes the Hubbard dimer with vanishing fundamental 
gap. 

The power functional inherits many of the problems of 
the Muller functional: There is no derivative discontinu¬ 
ity in the entire parameter range of the power functional 
except for the Hartree-Fock limit. In the low-interaction 
regime the solutions are weakly ferromagnetic. Like the 
Hartree-Fock approximation, the power functional ex¬ 
hibits an artificial abrupt magnetic transition with in¬ 
creasing interaction towards an antiferromagnetic con¬ 
figuration, albeit at a larger critical interaction. These 
states are intrinsically non-collinear. 

The absence of any derivative discontinuity also for 
insulating materials is expected to produce an artificial 


charge transfer between the constituents of large-scale 
electronic structure calculations. This cast severe doubt 
on the performance of such density-matrix functionals for 
complex systems. 

While the power functional lacks a derivative disconti¬ 
nuity, its chemical potential undergoes a continuous tran¬ 
sition between two linear functions, which has been ex¬ 
ploited to extract a band gap from data obtained further 
away from the integer particle number^®’^®^^®’®®. 

Our calculations indicate, however, that the band gap 
obtained from this extrapolation can be tuned by the free 
parameter a of the power functional between zero and 
the Hartree-Fock result. The band gap opens in non- 
collinear calculations only when in the antiferromagnetic 
regime, while it vanishes in the Miiller-type regime at low 
interactions. The opening of a band gap obtained by the 
extrapolation method and its tunability are features that 
persist in non-magnetic calculations, while the gap opens 
at a larger value of the power parameter than in the mag¬ 
netic calculation. These problems or signatures of them 
can be observed in previous calculations®^’^^’^®’®®’®®. 

The tunability of the band gap is similar to other meth¬ 
ods such as LDA-fU^ and hybrid density functionals®^. 
However, the latter methods exhibit a true derivative dis¬ 
continuity and their band gap does not shrink below the 
Kohn-Sham band gap, which is analogous to the non¬ 
interacting band gap of the Hubbard dimer. 

Approximations for ionization potentials®® and spec¬ 
tral functions®^ have been introduced on top of rDMFT. 
The latter method on the one hand yields spectra that 
agree well with experimental results for transition metal 
oxides®®’^®’®® for particular choices of the power func¬ 
tional parameter. On the other hand investigations on 
the Hubbard dimer®® suggest caution and claim that the 
underlying physics is not correctly treated. 

The problems presented here demonstrate potential 
fundamental flaws of the class of density-matrix func¬ 
tionals of this study. We hope that this study provides a 
useful reference point for the development of new density- 
matrix functionals. We believe furthermore that our 
findings call for new approaches for the construction of 
density-matrix functionals that make closer contact to 
the many-particle description of the electronic system^®. 
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